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We combine density functional theory (DFT) with molecular dynamics simulations based on an accurate 
atomistic force-field to calculate the pressure derivative of the melting temperature of magnesium oxide at 
ambient pressure - a quantity for which a serious disagreement between theory and experiment has existed for 
almost 15 years. We find reasonable agreement with previous DFT results and with a very recent experimental 
determination of the slope. We pay particular attention to areas of possible weakness in theoretical calculations 
and conclude that the longstanding discrepancy with experiment could only be explained by a dramatic failure 
of existing density functionals or by flaws in the original experiment. 



I. INTRODUCTION 

MgO is a major component of the Earth's mantle and 
so its thermodynamic properties at high pressures (P) 
and temperatures (T) are crucial to our understanding of 
its composition and evolution. It is arguably the simplest 
oxidei, being stable in the NaCl cubic structure at pres- 
sures up to at least 227 GPa at ambient temperature 2 -, 
and its simplicity and abundance make it a natural start- 
ing point for attempts to understand and model oxides 
of geophysical relevance^. 

Computer simulations based either on quantum me- 
chanics or on atomistic force-fields are playing an in- 
creasingly important role in geophysical research because 
the experimental difficulties at the extreme temperatures 
and pressures relevant to the Earth's mantle are consid- 
erable. However, for simulations to be of use, it is im- 
portant to be able to rely on the accuracy of theoretical 
descriptions of interactions between atoms and to under- 
stand their limitations. When bonding can be described 
accurately and efficiently, simulations allow many phys- 
ical properties of materials to be calculated at arbitrary 
temperatures and pressures. On the other hand, an in- 
ability to accurately calculate the properties of an oxide 
as simple as MgO would cast serious doubt on the suit- 
ability of computer simulations for quantitative studies of 
more complicated oxides such as (Mg,Fe)Si03-perovskite 
and (Mg,Fe)0 magnesiowiistite, which together make up 
about 90% of the lower mantle. 

Both MgO and MgSi03 are known to melt at temper- 
atures substantially above the geothermal profile, how- 
ever, a quantitative determination of their melting tem- 
peratures (T m ) at high pressures is a crucial parame- 
ter in rheological, geodynamical, and chemical differen- 
tiation models of the lower mantle^. Viscosity mod- 
els, for example, scale with the "homologous" tempera- 
ture (T/T m ), T being the actual temperature along the 
geotherm 6 . Chemical differentiation in the early, par- 
tially molten state of the mantle must have occurred at 
temperatures above the MgSiOa/MgO solidus, which is 
in turn determined by T m of the end-members. 



Until recently, only one experimental measurement 
of the melting temperature of MgO at high pressure 
existed^ and this extrapolated to a rather low value of 
5000 K for the melting temperature of MgO at core- 
mantle boundary pressures (130 GPa). If correct, this 
would imply that viscosity in the lower mantle is dom- 
inated by atomic diffusion in MgO — and suggest that 
partial melting may be the cause of the seismic anoma- 
lies at the bottom of the mantle^. However, atomistic 
modeling has consistently yielded a much steeper in- 
crease (dT m /dP) of the melting temperature of MgO 
with pressure^ - — The theoretical estimates of dT m /dP 
range from 88 K/GPaii to 270 K/GPe£, while Zerr and 
Boehler found a value of 36 K/GPa 7 . 

Because the melting slope is related through the 
Clapeyron relation (dT m /dP = T m AV/AE) to funda- 
mental physical properties of the material such as the 
change in molar volume upon melting (AV) and the la- 
tent heat AE, if the results of Zerr and Boehler were 
correct, it would point to a dramatic failure of atomistic 
models. However, even at low pressures, there are consid- 
erable difficulties associated with experimental measure- 
ments of the MgO melting points. Ronchi and Sheindlin 
have reported a zero pressure melting point of 3250 ± 20 
K— which differs significantly from the value of 3040±100 
measured by Zerr and Boehler. Further doubt has been 
cast on Zerr and Boehler's measurements by a very recent 
experiment in which Zhang and Fei have extrapolated a 
value of dT m /dP = 221 K/GPa from measurements of 
melting of (Mg,Fe)0 solid solutions at high pressure^. 

Here we combine molecular dynamics simulations with 
density functional theory (DFT) to determine the melt- 
ing slope of MgO. There have been two previous calcula- 
tions of the melting slope that relied heavily on DFT - one 
by AlfeiS and one by Aguado and Madden^. Our calcu- 
lations are intended to complement these studies and to 
demonstrate that, relative to the large discrepancy be- 
tween the calculated melting slope and that measured 
by Zerr and Boehler, there is agreement between values 
of the melting slope calculated by different groups using 
DFT. The calculations are also in better agreement with 
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the work of Zhang and Fei. 

We pay particular attention to analysing possible 
sources of error in our calculations and in previous cal- 
culations. Our calculations indicate that the rate of in- 
crease of the melting temperature with pressure is be- 
tween three and five times steeper than reported exper- 
imentally by Zerr and Boehler. Furthermore, this dis- 
crepancy does not appear to be explainable by statistical 
uncertainties in our calculations or by differences in the 
description of interatomic forces between the model po- 
tentials that we use for efficient statistical sampling and 
DFT. We are forced to conclude that either the local and 
generalised-gradient approximations to DFT fail spectac- 
ularly for solid and/or liquid MgO, or there are problems 
with the experimental results of Zerr and Boehler. Such 
a dramatic failure of density functionals for a material 
as simple as MgO would be very surprising and, to our 
knowledge, unprecedented. Therefore, given that the ex- 
periment of Zhang and Fei casts doubt on the results 
of Zerr and Boehler, it seems likely that current models 
of viscosity for the Earth's mantle which rely on these 
results need to be revised. 



II. PREVIOUS CALCULATIONS 

Many of the early calculations of the melting slope 
were performed using empirical atomistic models. These 
are energy functions of the atomic coordinates that 
were parametrized using low-temperature experimental 
or DFT data for crystalline MgO. There are several po- 
tential problems with these calculations. One poten- 
tial problem is that many of these atomistic models 
do not adequately describe electronic effects, such as 
ionic polarization, that may have a significant impact 
on thermodynamics. This is of particular concern for 
pairwisc-additive force-fields which do not contain any 
phenomenological representation of the response of elec- 
trons on an ion to changes in the ion's environment. A 
second potential problem is that the data to which these 
atomistic models were fit does not relate directly to the 
relevant thermodynamic (P, T) conditions and it doesn't 
relate directly to the liquid state. This means that one 
can be less confident of the models' applicability under 
these conditions where disorder and changes in volume 
may alter the electronic structure. Finally, if the quan- 
tity of data to which a model is fit is small it is relatively 
easy to achieve a good fit. However, one can never be 
sure that this fit results in a good underlying description 
of the forces on the ions. For these reasons, it is diffi- 
cult to assess the reliability of calculations of the melting 
slope that are based on these purely-empirical atomistic 
models. 

Parameter-free (or first-principles) approaches based 
on density functional theory and on the full description 
of the quantum electronic ground state have proven to be 
much more accurate and reliable than conventional force- 
fields for the calculation of the static and vibrational 



properties of crystalline MgO at low temperature^. A 
serious drawback of first-principles approaches, however, 
is their computational expense, which limits simulations 
to short time and length scales. Therefore, statistical 
sampling is usually poor and the precision with which 
thermodynamic properties can be calculated is low. Nev- 
ertheless, recent methodological advances and increasing 
computational resources have allowed the study of high- 
T thermodynamic properties of minerals in a few cases, 
including mcltingii and thermoelasticityi&. To find the 
reason for the discrepancy between theory and experi- 
ment on the melting slope of MgO, we will attempt to rule 
out as many of the possible reasons for this discrepancy as 
we can. Because simulations that rely solely on empirical 
or semi-empirical atomistic force-fields yield calculated 
melting slopes that differ by up to a factor of three and 
because their accuracy is very difficult to assess, we must 
assume, for the sake of the present argument, that they 
are untrustworthy. Therefore, we consider only the more 
recent calculations of dT m /dP that have been performed 
with substantial help from first- principles calculations. 

Alfe has calculated the melting slope of MgO with- 
out any reliance on atomistic force-fields by performing 
first principles molecular dynamics. He found a value 
for the melting slope of f02 ± 5 K/GPa 12 . However, we 
cannot rule out the possibility that his results are af- 
fected by uncertainties arising from short equilibration 
times and production runs. Aguado and Madden, on 
the other hand, have substantially reduced the proba- 
bility of poor equilibration and substantially increased 
the precision with which thermodynamic properties are 
calculated by using a highly-accurate atomistic potential 
that has been parametrized using DFT—. Equilibration 
and statistical sampling have been performed with this 
relatively-efficient force-field and they have used DFT to 
check the accuracy of the calculated energy differences 
between solid and liquid. They find a melting slope of 
125 K/GPa. Although the force-field that they use is 
very good, the configurations that they generate to cal- 
culate energy and volume differences can not be trusted 
as much as those that would be calculated if dynamics 
had been performed on the DFT potential energy surface. 
For example, they have parametrized their potential by 
fitting to DFT calculations of configurations from a num- 
ber of different solid phases. There is no guarantee that 
their force-field would be transferable to the liquid if the 
liquid structure differed strongly from these crystals. 

In this work, we calculate the melting slope using a 
similar approach to that of Aguado and Madden. We 
minimise finite-size effects and maximise the lengths of 
equilibration and production simulations by performing 
molecular dynamics with a highly-accurate and sophisti- 
cated atomistic force-field^. We use perturbation theory 
to correct the small differences between our force-field's 
description of the potential energy surface and that of 
DFT. We also take precautions to ensure that the con- 
figurations generated by our force-field are very close to 
those that would be generated directly from the DFT 
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potential energy surface: we parametrize this model by 
fitting to DFT forces, stresses, and energies calculated on 
the hot solid and the liquid; we also perform a first prin- 
ciples molecular dynamics simulation of liquid MgO to 
verify that we get a liquid structure that is very similar 
to that produced by our force-field. 

An important conclusion of the present work is that, on 
the scale of the discrepancy between theory and exper- 
iment, there is relative agreement between calculations 
of the DFT melting slope. It is important to note that 
there have been DFT-based calculations whose results 
differ strongly from those of Alfe, Aguado and Madden, 
and us. However, we draw a distinction between DFT- 
based calculations that simply use models that rely to 
some extent on DFT, and calculations of the DFT melt- 
ing slope. The latter are calculations which arrive at a 
close approximation of the melting slope that would be 
obtained from large scale molecular dynamics simulations 
on the DFT potential energy surface. As an example of 
the former kind of calculation we consider the study by 
Strachan et at who calculated the melting curve using a 
model that was fit to the DFT equations of state (equilib- 
rium volumes, bulk moduli, etc..) of the Bl and the high 
pressure B2 phases of MgO as well as the shear stresses 
along the transformation path between them. Therefore, 
this fit was to a very small amount of zero temperature 
DFT data and the resulting force-field was used directly 
to compute the melting line. 

The fit of our models is to high temperature solid and 
liquid DFT data and we converge this fit with respect to 
the quantity of DFT data (~ 5, 000 numbers are typically 
required). The closeness of the fit^ and the effectively- 
infinite amount of data used in the fit allows us to be 
confident that the force-field provides a very accurate 
description of the forces on the ions. However, the crucial 
point, as explained below, is that the role of our model is 
only to provide us with realistic statistically-independent 
hot solid and liquid configurations. The melting slope is 
computed by performing DFT calculations directly on 
these configurations so that, by first-order perturbation 
theory, we arrive at a close approximation to the DFT 
melting slope. Furthermore, we estimate the closeness of 
this approximation below. So, despite the fact that we 
have a very accurate DFT-parameterised model, our goal 
is not to calculate this model's melting slope but to use it 
as a stepping stone to calculate the DFT melting slope. 

Semi-empirical DFT calculations were performed by 
Cohen and Gong 8 , however, their Potential-Induced- 
Breathing (PIB) model imposes unphysical constraints 
on the density. For example, it is known that oxygen 
ions are highly polarisable but, within their approach, 
they remain spherically-symmetric. This results, among 
other effects, in a vast over-estimation of longitudinal op- 
tical phonon frequencies^. It is not known how oxygen 
polarisation affects the structure of liquid MgO, for ex- 
ample, but it is clear both from our classical and from 
our ab initio molecular dynamics simulations that oxy- 
gen ions acquire large dipoles in the disordered solid and 



in the liquid. 



III. CALCULATIONS 

We determine the melting slope dT m /dP of MgO at 
zero pressure by using the Clapeyron relation. We calcu- 
late Ay and AE with classical molecular dynamics and 
apply corrections to them using DFT. 

DFT calculations were carried out within the local den- 
sity approximation (LDA) using norm-conserving pseu- 
dopotentials with and without core corrections^ for Mg 
and O, respectively, and a plane wave basis set with 120 
Ry energy cut off. Simulation cells contained 64 atoms 
and the Brillouin zone was sampled with the T-point. 
Tests with 8 k-points yielded negligible (< 1%) differ- 
ences on solid-liquid energy differences, with respect to 
T-point sampling. 

In spite of their lower accuracy, model potentials can 
speed up considerably the task of calculating AV and 
AE from first principles if they are used as a "reference" 
model for the first-principles potential—. The model po- 
tential is used to generate statistically significant atomic 
configurations at the P-T conditions of interest and the 
first-principles values of AV and AE are then obtained 
by performing DFT calculations on those configurations 
only. We will show that, thanks to the quality of the 
model potential used in this work, the errors introduced 
by this procedure are significantly smaller than those in- 
trinsic to the standard approximations to DFT, which 
therefore remains the main source of uncertainty in our 
calculations. In order to achieve such a level of precision 
we use a model potential for MgO recently developed by 
us, which accounts for arbitrary aspherical distortions of 
the oxygen valence shells. Its parameters are obtained 
by best fit to DFT forces, stresses and energy in atomic 
configurations which are representative of the physical 
conditions of interest—. For this study we have used one 
potential ($j) which was optimized in the liquid at 3000 
K and P = GPa, and another (Q s ) which was opti- 
mized at the same P-T conditions in the solid. Average 
energies were set to be identical to the DFT values (this 
can trivially be imposed through an arbitrary additive 
constant). For both potentials, phonons, thermal expan- 
sion, and equations of state across a wide range of tem- 
peratures and pressures, are in very good agreement with 
experiments and with independent DFT calculations^. 

Fig. 1 shows AV and AE as extracted from long 
(~ 100 ps) molecular dynamics simulations of the solid 
(with $ s ) and the liquid (with $;) in a range of tempera- 
tures close to the experimental values for T m (3040 ± 100 
K£ or 3250 ± 20 K— ). Simulations were performed in 
cells containing 512 atoms under periodic boundary con- 
ditions. We verified that finite size effects on volume, 
compressibility, and thermal expansivity were negligible 
with this cell size. 

The first-principles values of AV and AE can be ob- 
tained by series expansion in the difference between the 
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reference potential and the first-principles potential 22 
The first-principles energy can be obtained as 

E{p = ($fp)fp = ($mp)mp + ($fp - $mp)mp + P(($fp ~ 

(^f P )m P )(*fp-$m P ))mp + 0(/3 2 ), where $ mp is the model 
potential ($ s or <!>;), <f>f p the first-principles potential, 
j3 = l/ksT (ks is the Boltzmann constant), and sta- 
tistical averages on the model or first-principles poten- 
tial are indicated with (...) mp and (...)fp, respectively. 
Similarly, for the first-principles value of the molar vol- 
ume we have, to lowest order in (3 and AV: Vt p ~ 
Ki P + V mp K T (P mp - Pf p } mp , where Pf p , mp is the pres- 
sure calculated from first-principles or with the model 
potential, and Kt is the isothermal compressibility. The 
first-principles value of AE (AV) is then obtained as the 
difference between the values of Ef p (Vf p ) in the liquid 
and in the solid. 

We computed E fp and V ip at 3070 K and P = GPa 
both in the liquid and in the solid by using twenty statis- 
tically independent configurations extracted from a long 
molecular dynamics run with the model potential. Each 
configuration was separated from the previous one by 
tens of picoseconds. Because the potentials have been ar- 
bitrarily given an energy offset so that (<J>f p — $ mp ) mp = 
at 3000 K, the first significant term of the series expan- 
sion for the energy is the linear term in f3. We verified 
that this term is indeed very small (-4.3 K in the liquid 
and -2.6 K in the solid), which implies that higher terms 
can be safely neglected. 

The same holds true for the volume, where we find 



that (P fp - P, 



mp/mp 



= 0.06 GPa, with a mean square 
deviation of 0.4 GPa, which means that uncertainties in 
the determination of the first-principles volumes are of 
the order of 1% . 

The very good performance of our model potential on 
the thermal expansion for the solid phased suggests that 
the agreement on AV between DFT and the potentials in 
Fig. 1 can be extended to all temperatures in the vicin- 
ity of 3070 K. A similar conclusion can be reached for 
the energy difference, based on the fact that energy fluc- 
tuations, and therefore heat capacities^, are correct to 
within 10%i£. We can summarize the above considera- 
tions by saying that the data of Fig. 1 represent the first- 
principles values of AE within f 0% and of AV within 2 
%. 

What is clear from Fig. f is that neither AV or AE 
is strongly temperature dependent and that the melting 
slope, therefore, depends approximately linearly on the 
melting temperature. The uncertainty in T m (and there- 
fore dT m /dP) is of the order of ~ 10% which is much 
less than the discrepancy with experiment on dT m /dP 
(> 300%) which we want to address in this work. Fits of 
straight lines to the data of Fig. 1 yield 



AE = 0.0295 
AV = 9.71 



1.97 x 10~ 7 T 
9.16 x 10~ 3 T 



(1) 



From this we can calculate melting slopes ranging from 
dT m /dP = 130 K/GPa if T m = 3050K to dT m /dP = 145 
if T m = 3250 K 



190 - 

O 180 -" 
5 

"5 170 - 
o> 

E 160 - 
o 

> 150 - 

140 — 



Liquid 
Solid 



-0.76 



-0.79 
2900 



- 


1 ' 












• Liquid 






■ Solid 










,1,1,- 



3000 3100 3200 

Temperature [Kelvin] 



3300 



FIG. 1. Energy and volume of liquid and solid MgO as a func- 
tion of temperature from molecular dynamics simulations. 
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FIG. 2. Pair correlation functions of solid and liquid MgO at 
~ 3100 K from molecular dynamics using a model potential 
compared to the results of first principles molecular dynam- 
ics of the liquid. Both our Car-Parrinello molecular dynamics 
(CPMD) simulations and Born-Oppenheimer molecular dy- 
namics simulations by Karki et al^ are presented. 



For the sake of completeness, we have attempted a 
determination of T m . This was achieved by first calcu- 
lating the melting temperature of by means of the 
two-phase method 9 , and then correcting the results us- 
ing perturbation theory^ 2 -. A simulation cell containing 
1024 atoms was used for the two-phase method. Pre- 
vious investigations*^ have concluded that the finite 
size effects are negligible with this size. We find that 
T m = 3010 ± 50 K for $/, the error reflecting only sta- 
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tistical uncertainties related to the two-phase method, 
and not other systematic errors due to the optimization 
procedure or to approximations to DFT. However, we 
caution that, in this simulation, potential is used to 
describe both the liquid and the solid phases, at variance 
with the calculations of AV^ and AE where two differ- 
ent potentials were employed to describe the liquid and 
the solid and which therefore resulted in a much more 
accurate calculation. From[T] therefore, we find that the 
value of the melting slope calculated with our model po- 
tential, before DFT corrections have been applied, is 127 
K/GPa. We can calculate DFT corrections to the Gibbs 
free energy difference between solid and liquid using the 
method described in Ref. HH. If we do so, we find that 
the DFT Gibbs free energy differences are 17% larger 
than those calculated with the model potential, indicat- 
ing that the model potential overstabilizes the liquid with 
respect to the DFT potential. This brings about a sim- 
ilar correction for T„^, so that the "DFT" value of T m 
is estimated to be about 3500 K, and the melting slope 
156 K/GPa. 

We have used DFT to correct the values of AV and 
AE obtained with our model potential and, as a result, 
can be confident that we have calculated the free energy 
difference between the solid and the liquid generated from 
$i with very close to DFT accuracy. However, there re- 
mains the possibility that the structure of the liquid gen- 
erated by is not realistic. Given the accuracy of this 
potential and that it has been parameterized from DFT 
calculations of the liquid, this possibility would appear to 
require a liquid-liquid phase transformation that involves 
changes in the electronic structure of the ions that cannot 
be captured by our model's phcnomcnological represen- 
tation of the electrons. It is important to consider this 
possibility, therefore, we have performed Car-Parrinello 
molecular dynamics (CPMD) simulations 2 ^ of 64 atoms 
of liquid MgO at 3050 K. These simulations were per- 
formed at zero pressure using variable cell dynamics 2 ^. 
A small fictitious mass of /i = 100 a.u. was used and 
temperature and pressure were corrected as described in 
Refl28l. The CPMD simulation was a continuation of 
a very long simulation using potential $/. After 2 ps 
of equilibration with CPMD, a 1.5 ps production run 
was used to compute the radial distribution functions of 
the liquid. The results are plotted in Fig[2] and com- 
pared with the results of Born-Oppenheimer molecular 
dynamics simulations carried out by Karki, Bhattarai, 
and Stixrude^. There is near-perfect agreement on the 
pair-distribution functions between the two independent 
first-principles simulations and our simulations of a 512- 
atom supercell using our model potential, $/. The agree- 
ment between the first principles simulations seems to 
rule out major structural artefacts of the starting con- 
figuration in both simulations. The agreement with the 
structure of the liquid generated with $/ confirms the 
reliability of this potential and also appears to rule out 
large finite-size effects on the pair-distribution functions. 
The results presented in FigJ5]are strong evidence that re- 



alistic liquid structures have been used to calculate DFT 
corrections to AV and AE and, therefore, that we have 
calculated these quantities with very close to DFT accu- 
racy. 

It is important to note at this point that our reported 
DFT results have been obtained within the LDA, as this 
approximation has proven to be very accurate in describ- 
ing low temperature properties. However, non local cor- 
rections to the LDA, such as those contained in general- 
ized gradient approximation theories (GGA), are known 
to have a significant effect on melting temperature o 12 i 29 
Therefore, as a test of the importance of exchange and 
correlation effects we repeated the analysis of energy dif- 
ferences with a GGA functional 30 . We find that average 
GGA energy fluctuations at 3000 K are within 12% of 
those calculated with the model potential. Moreover, we 
find that the correction to the Gibbs liquid-solid free en- 
ergy difference is only 2.7%, which implies a value for 
the GGA T m of 3090 K. This improves dramatically the 
agreement of T m with experiment, with respect to the 
LDA, and confirms that exchange and correlation effects 
are indeed important in the determination of T m . We 
caution, however, that the potential was constructed by 
fitting LDA quantities, so it is possible that the atomic 
configurations chosen for the comparison were not fully 
representative of the GGA potential. We did not attempt 
a determination of AV^ with GGA as this would require 
a very expensive equilibration with the GGA functional. 
AV would have to be an order of magnitude larger than 
the LDA value to resolve the discrepancy with experi- 
ment on the melting slope, which is highly unlikely. 

Alfe has pointed ou t 12 ' 31 that, because the energy gap 
between occupied and unoccupied electronic states is sig- 
nificantly smaller in the liquid than in the solid, it is more 
appropriate to perform DFT calculations with a finite 
electronic temperature. This results in a lowering of T m 
by approximately 500K. Consideration of this correction 
brings our LDA and GGA values of T m into very good 
agreement with those of Alfe. Alfe reports that the cor- 
rection to AE from this effect is almost 0.0036 a.u. per 
molecular unit. Because energy gaps calculated within 
LDA are generally too small (by as much as ~ 50%), the 
error incurred in AE is likely to be significantly smaller 
than this. In any case, for a fixed value of T m , the reduc- 
tion of AE by inclusion of the electronic entropy contri- 
bution to the free energy should increase the calculated 
melting slope thereby bringing the calculated value even 
further from the experimental value of Zerr and Boehler. 



IV. CONCLUSIONS 

To summarize, we find that that the DFT/LDA value 
for the melting slope of MgO ranges from ~ 130 K/GPa 
to ~ 150 K/GPa, depending primarily on the value cho- 
sen for T m and with an overall uncertainty of about 10- 
15% due to the model potential and to statistical sam- 
pling. We can safely conclude that the DFT/LDA result 
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is a factor of 3 to 4 larger than the value of 36 K/ GPa 
found in Zcrr and Boehler's experiment. 

There remains a small difference in the melting slopes 
calculated here and by Alfe. The source of this difference 
is unknown, but may be due to differences in the details 
of the DFT calculations such as our different pseudopo- 
tential representations of the inert core electrons. The 
important point is that this difference is substantially 
smaller than the discrepancy with the experimental re- 
sults of Zerr and Bochlcr. The present work, when taken 
together with the results of Alfe and Aguado and Mad- 
den, shows that on the scale of the discrepancy between 
calculations and experiment, there appears now to be 
a convergence of the results of calculations of the melt- 
ing slope. The only major source of error that would 
effect all of these calculations is the approximation to 
the exchange-correlation energy. However, tests with 
a GGA exchange-correlation functional show that, al- 
though the choice of functional changes the values of T m 
and dT m /dP, these changes are relatively small. There- 
fore it seems highly unlikely that inadequacies of the DFT 
approximations used can fully explain the historical dis- 
crepancy between theory and experiment. 

The possibility of problems with the experiment of 
Zerr and Boehler— have previously been suggested^^ 
and the fact that there is a large disagreement be- 
tween experimental measurements of both the melting 
temperature 7 -^ and the melting slope 7 -^ suggests that 
more experimental work is necessary. A determination of 
the density change and latent heat at zero pressure are 
crucial to resolve the issue. The disagreement between 
our DFT results and experiment could also be explained 
by assuming that the slope of the melting curve is ini- 
tially very steep, but that it flattens out very quickly, 
perhaps due to a liquid structure which changes rapidly 
under pressure to being much more similar to the solid. 
However, this explanation would not be compatible with 
the findings of Alfe, or Aguado and Madden, who explic- 
itly computed the melting temperature at high pressures. 
Therefore, we must discount this possibility. 

We note that the suggestion by Aguado and Madden 
that the discrepancy between theory and experiment can 
be explained by the existence of a solid phase with a 
lower free energy than the rocksalt structure cannot be 
correct. This is because, among possible solid phases, the 
one with the lowest Gibbs free energy has the highest 
melting temperature. Therefore, the melting curve for 
a more stable solid phase should lie above the rocksalt 
melting curve, i.e. it should have a higher T m at every 
pressure P. 

The geophysical implications of a steep melting slope 
for MgO are manyfol d 12 ' 15 . A steep melting slope implies 
that the melting temperature of (Mg,Fe)0, of which 
MgO is an end-member, is likely to be substantially 
higher than the geotherm and comparable to the melting 
temperature of (Mg,Fe)Si03 at lower mantle conditions. 
This suggests, among other consequences, that large 
scale melting may never have occurred in the mantle 5 . 
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